#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <assert.h>

#define N 100
#define EPS 1E-6
#define PI 3.141592654

/********** 2d geometry utils **********/

typedef struct {
  double x, y;
} point;

typedef struct {
  point p, q;
} line;

static inline double frand(double a, double b) {
  return a + (b - a) * rand() / RAND_MAX;
}

static inline int sign(double x) {
  return x < -EPS ? -1 : x > EPS;
}

static double cross(point a, point b, point c) {
  return (b.x - a.x) * (c.y - b.y) - (b.y - a.y) * (c.x - b.x);
}

static int intersects(point a, point b, point c, point d) {
  return sign(cross(a, b, c)) * sign(cross(a, b, d)) < 0 &&
         sign(cross(c, d, a)) * sign(cross(c, d, b)) < 0;
}

static double dist(point a, point b) {
  return hypot(a.x - b.x, a.y - b.y);
}

static double max(double x, double y) {
  return x > y ? x : y;
}

static double min(double x, double y) {
  return x < y ? x : y;
}

/*
 * generates a planar graph
 * stores into adj matrix g[N][N]
 */
void gen(int (*g)[N], int n, int seed) {
  static point p[N];
  static int deg[N];
  int i, j, x, y, ok, fails = 0;

  if (seed != -1)
    srand(seed);

  for (i = 0; i < n; i++) {
    p[i].x = frand(0, 1);
    p[i].y = frand(0, 1);
    deg[i] = 0;
    for (j = 0; j < n; j++)
      g[i][j] = 0;
  }

  while (fails < n * n) {
    do {
      x = rand() % n;
      y = rand() % n;
    } while (x == y || g[x][y] || (deg[x] > 3 && deg[y] > 3) || dist(p[x], p[y]) > 2.0 / sqrt(n));

    ok = 1;
    for (i = 0; ok && i < n; i++)
      for (j = 0; ok && j < i; j++) if (g[i][j])
        if (intersects(p[x], p[y], p[i], p[j]))
          ok = 0;
    if (ok) {
      g[x][y] = g[y][x] = 1;
      deg[x]++;
      deg[y]++;
      fails = 0;
    } else
      fails++;
  }
  puts("gen completed.");
}

/* check if it is a planar graph
 * return number of intersecting edge pairs 
 */
int check(int (*g)[N], int n, const point *p, int verbose) {
  static const point *edge[N * 3][2];
  int x, y, i, j, iss = 0, m = 0;
  for (i = 0; i < n; i++)
    for (j = 0; j < i; j++) if (g[i][j]) {
      edge[m][0] = &p[i];
      edge[m++][1] = &p[j];
    }
  for (x = 0; x < m; x++)
    for (y = 0; y < x; y++)
      if (intersects(*edge[x][0], *edge[x][1], *edge[y][0], *edge[y][1]))
        iss++;
  if (verbose) {
    if (iss == 0)
      puts("pass");
    else
      printf("intersected edge pairs: %d\n", iss);
  }
  return iss;
}

/*
 * randomly shuffles to points
 */
void place_random(point *p, int n, int seed) {
  int i;
  if (seed != -1)
    srand(seed);
  for (i = 0; i < n; i++) {
    p[i].x = frand(0, 1);
    p[i].y = frand(0, 1);
  }
}

/*
 * place the points circularly
 */
void place_circular(point *p, int n) {
  int i;
  for (i = 0; i < n; i++) {
    p[i].x = 0.5 + 0.5 * cos(i * PI * 2 / n);
    p[i].y = 0.5 + 0.5 * sin(i * PI * 2 / n);
  }
}

/*
 * normalizes the points
 * fit them into the unit box (0, 1) x (0, 1)
 */
void normalize(point *p, int n) {
  int i;
  double xmin, ymin, xmax, ymax, dx, dy;

  xmin = p[0].x;
  ymin = p[0].y;
  xmax = p[0].x;
  ymax = p[0].y;
  for (i = 0; i < n; i++) {
    if (xmin > p[i].x)
      xmin = p[i].x;
    if (ymin > p[i].y)
      ymin = p[i].y;
    if (xmax < p[i].x)
      xmax = p[i].x;
    if (ymax < p[i].y)
      ymax = p[i].y;
  }
  dx = xmax - xmin;
  dy = ymax - ymin;
  if (0 < xmin && xmax < 1 && 0 < ymin && ymin < 1 && dx > 0.5 && dy > 0.5)
    return;
  for (i = 0; i < n; i++) {
    p[i].x -= xmin;
    p[i].y -= ymin;
    p[i].x /= dx;
    p[i].y /= dy;
  }
  for (i = 0; i < n; i++) {
    assert(!isnan(p[i].x));
    assert(!isnan(p[i].y));
  }
}

static double max_speed(const point *v, int n) {
  double res, tmp;
  int i;
  for (i = 0; i < n; i++) {
    tmp = hypot(v[i].x, v[i].y);
    if (i == 0 || res < tmp)
      res = tmp;
  }
  return res;
}

/********** json **********/

void json_write(int (*g)[N], int n, const point *p, FILE *file, int id) {
  int i;
  /* points */
  fprintf(file, "    [\n");
  for (i = 0; i < n; i++)
    fprintf(file, "      {x: %lf, y: %lf},\n", p[i].x, p[i].y);
  fprintf(file, "    ],\n");
}

void json_write_begin(int (*g)[N], int n, const point *p, FILE *file) {
  int i, j;
  fprintf(file, "model = {\n");
  fprintf(file, "  n: %d,\n", n);
  fprintf(file, "  g: [\n");
  for (i = 0; i < n; i++) {
    for (j = 0; j < i; j++) if (g[i][j])
      fprintf(file, "    {i: %d, j: %d},\n", i, j);
  }
  fprintf(file, "  ],\n");
  fprintf(file, "  ps: [\n");
}

void json_write_end(int (*g)[N], int n, const point *p, FILE *file) {
  fprintf(file, "  ]\n");
  fprintf(file, "}\n");
}

/********** solver **********/

void solve(int (*g)[N], int n, point *p, int iter, const char *filename) {
  static point v[N];
  static point q[N];
  const double K2 = 0.0005;
  double K1 = 0.05;
  int i, j, id = 0;
  int error, rate1, rate2;
  double dt, vmax;
  FILE *file;

  puts("solving...");

  /* init point position */
  place_circular(p, n);

  /* write begin */
  file = fopen(filename, "w");
  json_write_begin(g, n, p, file);
  json_write(g, n, p, file, id++);

  /* iteration */
  for (i = 0; i < n; i++)
    v[i].x = v[i].y = 0;
  rate1 = rate2 = n * n;
  while (iter--) {
    for (i = 0; i < n; i++) {
      for (j = 0; j < n; j++) {
        dt = dist(p[i], p[j]);
        dt = dt * dt;
        if (g[i][j]) {
          v[i].x += (p[j].x - p[i].x) * dt * K1;
          v[i].y += (p[j].y - p[i].y) * dt * K1;
        } else {
          dt = max(dt, 0.01);
          v[i].x -= (p[j].x - p[i].x) / dt * K2;
          v[i].y -= (p[j].y - p[i].y) / dt * K2;
        }
      }
    }
    for (i = 0; i < n; i++) {
      q[i].x = p[i].x;
      q[i].y = p[i].y;
      p[i].x += v[i].x;
      p[i].y += v[i].y;
    }
    normalize(p, n);
    for (i = 0; i < n; i++) {
      v[i].x = p[i].x - q[i].x;
      v[i].y = p[i].y - q[i].y;
    }

    /* write json */
    if (iter % 5 == 0)
      json_write(g, n, p, file, id++);

    /* statistics */
    error = check(g, n, p, 0);
    rate1 = (rate1 * 50 + error * 50) / 100;
    rate2 = (rate2 * 10 + error * 90) / 100;
    vmax = max_speed(v, n);
    printf("iter %d  rate1 = %d, rate2 = %d, vmax = %lf, K1 = %lf\033[K\r", iter, rate1, rate2, vmax, K1);
    fflush(stdout);
    if (vmax < 0.005 && abs(rate1 - rate2) <= 5) {
      if (vmax < 0.0001)
        break;
      else
        K1 = min(4.0, min(K1 * 1.5, K1 + 0.2));
    }
  }
  puts("");

  /* write end */
  json_write_end(g, n, p, file);
  fclose(file);
  printf("written to `%s'\n", filename);

  puts("solve completed.");
}

/********** main **********/

int main() {
  static int g[N][N];
  static point p[N];
  int n, seed;
  gen(g, n = N, seed = time(0));
  solve(g, n, p, 1000, "draw_data.js");
  check(g, n, p, 1);
  return 0;
}
